# -*- coding: utf-8 -*-
# @Author  : longbhu
# @Time    : 2025/3/12 13:38
# @Function: 双线性的插值方法
"""

"""
import os
import rasterio
from rasterio.warp import reproject, Resampling
import numpy as np

from calc_v1.tif_tools import tif_mask_by_shapefile


# from pathlib import Path

def interpolate_to_target_resolution(input_path: str, output_path: str = None, target_resolution: float = 0.01) -> None:
    """
    将输入的 TIF 文件插值到目标分辨率（默认 0.01°），并保存结果。

    参数:
        input_path (str): 输入的 TIF 文件路径。
        output_path (str): 输出的 TIF 文件路径。
        target_resolution (float): 目标分辨率（单位：度），默认为 0.01°。

    返回:
        None
    """
    # 自动生成输出路径
    if output_path is None:
        base_name = os.path.splitext(input_path)[0]
        ext = os.path.splitext(input_path)[1]
        meters = int(target_resolution * 111000)  # 1度≈111km
        output_path = f"{base_name}_{meters}m{ext}"

    # 打开输入文件
    with rasterio.open(input_path) as src:
        # 读取数据
        data = src.read(1)
        transform = src.transform
        crs = src.crs
        nodata = src.nodata  # 获取无效值

        # 判断 CRS 是否缺失
        if crs is None:
            print("警告：输入文件未包含有效的 CRS，将使用默认 CRS (EPSG:4326/WGS84)。")
            crs = "EPSG:4326"  # 设置默认 CRS 为 WGS84

        # 如果 nodata 值不是 -9999，将其修改为 -9999
        if nodata is not None and nodata != -9999:
            data[data == nodata] = -9999  # 替换 nodata 值为 -9999
            nodata = -9999.0  # 更新 nodata 值

        # 将无效值替换为 nodata
        data = data.astype('float32')  # 转换为浮点型以支持 NaN
        nodata = float(nodata) if nodata is not None else -9999.0  # 统一为浮点数
        # 使用浮点数阈值判断（1e-6为误差容忍度）
        mask = np.isclose(data, nodata, atol=1e-6, rtol=0)
        data[mask] = np.nan


        # 计算目标尺寸
        target_width = int((src.bounds.right - src.bounds.left) / target_resolution)
        target_height = int((src.bounds.top - src.bounds.bottom) / target_resolution)

        # 创建目标数组
        target_data = np.full((target_height, target_width), nodata, dtype=data.dtype)
        print(f"创建目标数组完毕！{target_data.shape}")

        # 计算目标变换矩阵
        target_transform = rasterio.Affine(
            target_resolution, 0, src.bounds.left,
            0, -target_resolution, src.bounds.top
        )

        # 使用双线性插值进行重投影
        reproject(
            source=data,
            destination=target_data,
            src_transform=transform,
            src_crs=crs,
            dst_transform=target_transform,
            dst_crs=crs,
            resampling=Resampling.bilinear,  # 使用双线性插值
            src_nodata=nodata,  # 指定源数据的无效值
            dst_nodata=nodata  # 指定目标数据的无效值
        )

        # 将 NaN 替换回原始无效值（如果存在）
        target_data[np.isnan(target_data)] = nodata
        print("将 NaN 替换回原始无效值（如果存在）")

        # 保存插值结果
        with rasterio.open(
            output_path,
            'w',
            driver='GTiff',
            height=target_height,
            width=target_width,
            count=1,
            dtype=target_data.dtype,
            crs=crs,
            transform=target_transform,
            nodata=nodata,  # 设置输出文件的无效值
            # compress='lzw',  # 压缩格式
            # bigtiff='IF_SAFER'  # 防止内存溢出
        ) as dst:
            dst.write(target_data, 1)


            # 显式生成统计信息（关键修复）
            # dst.write_colormap(1, {})  # 清空伪彩色表
            dst.update_tags(ns='rio_stats', STATISTICS_VALID_PERCENT=100)  # 标记有效统计
            # 强制重新计算统计信息（排除无效值）
            # stats = dst.stats(1)
            stats = dst.statistics(1)  # 或者使用正确的属性名
            print(f"实际数据范围: {stats.min} ~ {stats.max}")  # 验证输出范围

    print(f"Interpolated data saved to {output_path}")

# 示例调用
# input_path = r'F:\test_files\插值验证程序的输入数据\双线性插值\2023年降水量空间分布栅格图.tif'
# output_path = r'F:\test_files\插值验证程序的输入数据\双线性插值\2023年降水量空间分布栅格图_30m.tif'
# # 插值到 0.00027777778° 分辨率
# interpolate_to_target_resolution(input_path, output_path, target_resolution=0.0027777778)


# input_path = r'G:\GEP_data\wr\全年降雨量\2023_precip_china.tif'
# # output_path = r'F:\test_files\插值验证程序的输入数据\双线性插值\2023年降水量空间分布栅格图_30m.tif'
# # 插值到 0.00027777778° 分辨率
# interpolate_to_target_resolution(input_path, target_resolution=0.00027777778)

input_folder = r'G:\需要插值的变量'
for file in os.listdir(input_folder):
    if file.endswith('.tif'):
        input_path = os.path.join(input_folder, file)
        output_folder = os.path.join(input_folder, '30m')
        os.makedirs(output_folder, exist_ok=True)
        output_path = os.path.join(output_folder, file.replace('.tif', '_30m.tif'))

        interpolate_to_target_resolution(input_path, output_path, target_resolution=0.007)
        print(f"Interpolated data saved to {output_path}")

shp_input_folder = os.path.join(input_folder, '30m')
shapefile_folder = r'E:\shp\2025购买边界shp2021\省级'
for file in os.listdir(shp_input_folder):
    if file.endswith('.tif'):
        input_path = os.path.join(shp_input_folder, file)
        new_folder = file.split('_')[0]
        output_folder = os.path.join(shp_input_folder, new_folder)
        os.makedirs(output_folder, exist_ok=True)

        tif_mask_by_shapefile(input_path, shapefile_folder, output_folder, city_name='xzqdm')  # xiang_name
        print(f"{input_path}, Shp masked data saved to {output_folder}")






